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Abstract. High harmonic generation in polarizable multi-electron systems is 
investigated in the framework of multi-configuration time-dependent Hartree-Fock. 
The harmonic spectra exhibit two cut offs. The first cut off is in agreement with the well 
established, single active electron cut off law. The second cut off presents a signature of 
multi-electron dynamics. The strong laser field excites non-linear plasmon oscillations. 
Electrons that are ionized from one of the multi-plasmon states and recombine to the 
ground state gain additional energy, thereby creating the second plateau. 
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When an intense laser pulse is focused onto a noble gas jet, high harmonic 
generation (HHG) takes place. High harmonic radiation is created in a three step 
process pQ. The valence electron is set free by tunnel ionization. In the continuum, 
the electron is accelerated and follows the quiver motion of the laser field. When the 
laser field changes sign, the electron is driven back towards the parent ion. Finally, the 
electron recombines to the ground state upon recollision, and an xuv photon is emitted. 

The theory of HHG is based on the single-active-electron (SAE) approximation [2|, 
assuming that only the valence electron interacts with the strong laser field while the 
residual electron core remains frozen. The valence electron and the core electrons are 
regarded as uncorrelated. HHG has been performed with noble gas atoms and clusters 
[3], and with small molecules [UEJ- All experiments were found to be in agreement with 
SAE theory. 

Experimental [El E] and theoretical [Oj evidence was found that SAE theories 
cannot describe optical field ionization of highly polarizable systems, such as large 
molecules and metallic clusters. Due to the high electron mobility and polarizability, 
a factorization into valence and core electrons is no longer valid and the complete, 
correlated multi-electron (ME) dynamics has to be taken into account. This raises the 
question as to which extent the SAE approximation is applicable to non-perturbative 
phenomena in complex materials fUJ • 

In this article we investigate HHG in highly polarizable molecules by an 
one-dimensional (ID) multi-configuration time-dependent Hartree-Fock (MCTDHF) 
analysis. MCTDHF is a recently developed method allowing to account for the electron 
correlation in a numerically converged manner [T2*l [T3\ IT^ \W\ ITo*] . Our analysis reveals 
that in contrast to HHG in noble gases, where the harmonic spectrum exhibits one 
plateau and cut off, a second cut off is identified, extending far beyond the standard 
cut off. This second cut off originates from the ME nature of the bound electrons. The 
strong laser field excites non-linear, collective electron oscillations and populates multi- 
plasmon states that oscillate at a multiple of the plasmon frequency. The second plateau 
is generated by electrons that ionize from the multi-plasmon states and recombine to 
the ground state. The energy difference between excited and ground state determines 
the difference between first and second cut off. 



1. The MCTDHF method 



Here, we demonstrate the general idea of the MCTDHF-ansatz by means of an example, 
containing two ID particles. For simplicity we will not take spin into consideration, 
although it is included in the calculation presented below. For an extensive review of 
the MCTDHF theory and formalism we refer to [15. and references therein. 
MCTDHF makes the ansatz 

i m 

*(x,y,t) = —= A jlJa (t) [^p h {x-t)ip j2 {y-t) - ^ j2 (x;t) Vjl (y;t)] (1) 
Thus MCTDHF consists in approximating the exact wave function as linear combination 
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Figure 1. Approximation of a 2-dimensional, correlated wave function ^>(x,y;t) as 
a sum of 6 determinants, indicated as rectangular patches. The expansion functions 
tpi belonging to their corresponding determinant are drawn along the axes in their 
respective color. 



of (^j different Slater determinants. Since the coefficients and the m linearly 
independent expansion functions are time-dependent, an additional constraint is needed. 
Without loss of generality, it is most convenient to impose ortho-normality on the 
expansion functions, i.e. {<f%{t)\ifij{t)) = 8ij. Additionally, we impose (^■^fi(t)\'Pj(t)\ = 
to uniquely define the expansion in (0) [TK] . 

To derive equations of motion the Dirac-Frenkel variational principle f?l ITS] . 
Sty i-^z — H(t) ty \ — 0, is imposed. Here, H denotes the time-dependent Hamiltonian 



describing the electronic system. The such derived equations are of the form 

i±A = F[<p]A (2) 

i±- t p = G[A,<p)<p. (3) 

The time evolution is obtained by applying non-linear operators F and G on vectors 
A = (Ai t 2, Ai-i,n) T and y> = (ip%, (p m ) T , respectively. For details on the equations 
see [To] . 

The principle of MCTDHF can be seen in figure [0 The fully correlated wave 
function ty(x,y;t) is approximated by 6 Slater determinants, implying 4 expansion 
functions per particle. Thus the wave packet is reassembled in a kind of "patch work". 
Each patch represents one single Slater determinant. In figure ^ we have, for simplicity, 
assumed the single electron orbitals to be of rectangular shape. However, we emphasize 
that generally, the expansion orbitals will be a priori unknown. It is important to note 
that not only the expansion coefficients evolve with time, but the single particle orbital, 
too. Hence MCTDHF may be interpreted as a truncated configuration-interaction (CI) 
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expansion, in which both, coefficients, Aj 1 j 2 (t), and orbitals, (fj(x; t), are optimized. For 
every time step and fixed number of configuration an optimal expansion is warranted 
by the Dirac-Frenkel variational principle. Thus resulting in a compact representation 
of the wave function and hence compressing the necessary storage amount. 

MCTDHF is applicable to more complex systems because already with "small" 
configuration numbers the "essential physics" is covered, i.e. the number of physically 
important expansion orbitals, is always much smaller than the number of time- 
independent basis functions in conventional approaches. Therefore MCTDHF scales 
more slowly, allowing the treatment of small molecules beyond state of the art 2- 
electrons-calculation. 

By increasing the number of configurations MCTDHF allows a systematic inclusion 
of correlation, converging to the exact solution for m — > oo. 

These advantages are achieved at the expense of linearity and locality, since the 
evolution equations in MCTDHF are both, non-linear and non-local. Just as time- 
dependent Hartree-Fock (TDHF), also MCTDHF with a finite number of configurations 
suffers, in principle, from the problem of non-linearity of the evolution equations which 
may lead to a violation of the superposition principle. However, in the case of MCTDHF 
this problem is greatly reduced by adding additional configurations. 

The ability of MCTDHF to correctly describe the correlated dynamics of electrons 
has been assessed in recent studies jlHl EI| using two-electron systems. 

Here, we report calculations for n = 4 electrons using m = 8 expansion functions. 
The Schrodinger equation is solved in a simulation box with size / = ±360 at. u. and 
on a uniform ID-grid with 2400 grid points, using a second order finite difference 
representation. To avoid reflection at the boundaries, complex absorption potentials 
(CAP) are used. That is, the total Hamiltonian is modified by adding iX^O ~~ 
cos [ir{xi — x ca p)/(\l\ — a; cap )]} for \xi\ > x cap = 270at.u. and otherwise. The time- 
integration is performed by a self-adaptive, high-order Runge-Kutta integrator with 
a relative numerical accuracy of 10~ 8 . Convergency was checked with respect to all of 
these parameters. In particular, increasing the number of expansion functions to m = 12 
does not change our finding and changes for instance the ionization yield by less than 
4%, indicating that our calculations are essentially converged. 



2. High harmonic generation in polarizable molecules 

The ID MCTDHF analysis is based on the solution of the ID Schrodinger equation 
for the n = 4-electron potential V = J27=i ~ Vn{ x i) +Ej>iK(2 ; ! ~ x j) ■ Here, V n = 



Zj ' \Jx\ + a\ refers to the nuclear binding potential, Z is the charge state, and a n is the 
shielding parameter of the electron-nucleus interaction. The shielded model potential 
represents an atom or a small cluster. We believe that it is closer to a small cluster, 
with a harmonic oscillator potential part close to the center, and a Coulomb far-range 



potential far away from the center. Further, V e = 1/J (xi — Xj) 2 + a\ represents the 
electron-electron interaction potential with shielding parameter a e . The laser is coupled 
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Figure 2. Spectra of the dipole moment squared, | d(w) | 2 , of a highly polarizable 
(ao = 31 A 3 ), 4-electron model-system, I p = 0.5 at. u., (thick full line), and 
corresponding SAE-calculation for the same I p (thin dashed line). The lower x- 
axis gives the harmonic order, the upper x-axis gives the classical return energy, 
E T = (E—I p )/Up with E the harmonic photon energy. The standard cut off harmonic is 
at [E — Ip)/Up = 3.17 and is marked with an arrow. Laser parameters: Ao = 1000 nm, 
peak intensity I = 2 x 10 14 W/cm 2 , FWHM pulse duration r = 4T , optical period 
T = 3.33 fs, Gaussian envelope. 



in velocity gauge and in dipole approximation. Atomic units are used throughout, unless 
otherwise stated. 

ID ME simulations tend to overestimate the polarizability. To keep the 
polarizability at a reasonable level, the ionization potential had to be chosen slightly 
higher than usual values of complex materials (for instance benzene: J p = 0.35 at. u.). 
The softening parameter used for the SAE system is a n = 1.414, and the parameters 
to model highly polarizable atoms are a n = 0.80, a e = 1.0, and Z = 4. The binding 
energy of the 4-electron ground state is = 8.5 at.u. and the successive single electron 
ionization potentials are given by 0.5, 1.07, 3.09 and 3.93 at.u. The static polarizability 
is calculated by using the relation a = 1/8, J Ap(x)xdx, where Ap is the change in 
electron density caused by the field £. We find a polarizability of «o = 31 A 3 , which 
lies between the polarizability of transition metal atoms and clusters, for example, Nb: 
a = 15 A 3 , Ip = 0.248 at.u.; C 60 : a = 80 A 3 , J p = 0.279 at.u. Finally, the laser 
parameters are: center wavelength A = 1000 nm, peak intensity I = 2 x 10 14 W/cm 2 , 
Gaussian envelope with FWHM width r = 4T , and oscillation period T = 3.33 fs. The 
evolution of the wave function is calculated between 20 optical cycles before and 80 
optical cycles after the laser pulse maximum. 
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Figure 3. Contour plot of a time- frequency analysis of the 4-electron spectrum in 
Figure. [3 (right panel). The window-function is a Gaussian pulse with 0.2 optical 
cycles FWHM duration. Contours differ by a factor of 10 3 , decreasing from left to 
right. The return time t T is plotted versus the harmonic frequency normalized to the 
laser frequency. Here the thick dashed line corresponds to the standard cut off, while 
the thick full line represents the second cut off. The left panel shows the corresponding 
laser electric field. The time of birth and the time of return for a classical electron 
acquiring the maximum kinetic energy during its excursion in the laser field are marked 
by a dot and a cross, respectively. 

In figure |2] the harmonic spectrum is shown for a 4-electron system (full line) and 
a SAE system (dashed line) with the same HOMO (highest occupied molecular orbital) 
ionization potential I p = 0.5 at.u. The harmonic spectrum is obtained as the modulus 
of the Fourier transform of the dipole expectation value, d{t) = (^(t)\Y^i=iXi\^{t)). 
d(t) was sampled 256 times per optical cycle. Note, however, that the time step in- 
between these sample points was self-adaptive. The cut off energy E = I p + 3.17£/ p is in 
agreement with the standard SAE cut off law [HI2]- Here U p = (8, /2u ) 2 = 0.68 at.u. 
is the ponderomotive energy, £ is the laser peak field strength, and u denotes the 
laser circular center frequency. The ME spectrum reveals in addition to the regular, 
first cut off a second one. Note, that in the SAE system, ionization is saturated before 
the peak of the laser pulse, which is not the case for the multi-electron case, where 
the saturation intensity is increased due to the molecule's polarizability Ej. Thus 
the early saturation of the SAE system reduces the probability of electron trajectories 
that return with high energy and therewith results in a low high harmonic yield of the 
plateau. 

To identify the origin of the second plateau we have performed a time-frequency 
analysis of the ME spectrum, depicted in figure El The dipole moment is truncated 
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by the window function a w (t;T w ) = 1/(ttT w ) 1 / 4 exp [-{t - t r ) 2 /(2T*)} with T w = 0.2T 
and then Fourier transformed, i.e. d(ui,t T ;T w ) = ^{a w {t v \ T w )d(t T )}. The harmonics 
corresponding to the first and second plateau are depicted by the thick, dashed, and 
the thick, full lines, respectively. The time- frequency analysis is a way to connect the 
quantum mechanical result with the classical three-step model JTj model of HHG. It cuts 
small chunks in time out of the wavefunction and determines their energy at the time of 
return to the parent ion. For a SAE harmonic spectrum, the resulting graph of harmonic 
energy versus return time is very close to the result obtained by the classical three-step 
analysis. This correspondence allows an interpretation of the the time-frequency plot 
and of HHG in terms of classical trajectories. 

The thick dashed line in figure H3 denotes the regular, first cut off. A comparison 
to the SAE three-step model allows us to determine the importance of ME effects in 
HHG. The electron return phase of the cut off trajectory creating the highest harmonic 
in figure El is around 60° after the pulse maximum. This is shifted with respect to the 
three-step model cut off trajectory that is born at 163° before and returns at 80° after 
the pulse maximum pP (see left panel in figure EJ). The difference in the return times 
arises from a many-body effect. In ME systems the laser field induces a polarization of 
the bound electrons that exerts a repelling force on the ionizing electron. This additional 
potential decreases rapidly with the distance from the parent system and hence, affects 
the electron trajectories only in the vicinity of the nucleus. Therefore, it mostly shifts 
the birth and return time of the electron trajectories and only weakly affects the highest 
achievable cut off energy. 

The thick full line in figure El corresponds to the second cut off. Surprisingly, the 
first and second plateau show similar patterns. The maximum energy in each half cycle 
occurs at the same return phase for both plateaus. This strongly indicates that the 
harmonics in both plateaus are generated by the same electron trajectories, starting 
from different initial states. The strong laser field brings the medium in a coherent 
superposition of ground and excited states. HHG from excited states can take place as 
long as the phase of ground and excited states are coherently locked. 

In ME systems there exist two types of excitation, collective excitations and 
individual particle excitations. Single electron excitations can be excluded for the 
following reasons: (i) The SAE calculation in figure El does not show a second plateau, 
(ii) The energy difference between the first and second cut off is larger than the HOMO 
ionization potential, (iii) The absence of doubly ionized states excludes HHG from a 
deeper bound electron, (iv) While HHG from a coherent superposition of the ground 
state and an excited single-electron state does indeed produce two plateaus, it does not 
result in an over all increase of the standard cut off law [TUj . because the HOMO electron 
is born with the same energy after ionization regardless of its initial state. 

The excitation mechanism responsible for the second plateau is revealed in figure |U 
In sub-figure Ufa) the dipole moment, d(t), is plotted as a function of time. We find 
that the dipole moment exhibits oscillations even after the laser pulse, proving a laser 
induced excitation of the system. The close-up in the inset of figure Ufa) shows that the 
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Figure 4. (a) Time dependence of the applied laser electric field (lower panel) and 
the resulting dipole moment, d(t) (upper panel), of the highly polarizable molecule in 
figure |21 In the inset we show a magnified part of the dipole moment to illustrated the 
remaining excitations after the end of the laser pulse, (b) Maximum excursion of the 
center of gravity of the electron density (upper panel) and the phase shift, </>, between 
the dipole signal and the laser electric field (lower panel) as a function of the applied 
laser frequency. The full line is obtained for to = 8 while the dots correspond to 
to = 12 calculations. The closeness of agreement for both sets of calculations confirms 
that the electronic structure of the model molecule is well approximated and essentially 
converged. In (b) we have used a continues-wave laser which was linearly switched on 
over 4 optical cycles reaching a maximum intensity of 2 x 10 13 W/cm 2 . The pulse was 
propagated over 13 optical cycles. 



excitation oscillates at a frequency u p = 0.35 at.u. In figure^ b) the excitation spectrum 
of the system is determined by probing the response of the system to a plane wave laser 
signal as function of the laser center frequency uq. The maximum excursion of the 
center of gravity of the electron density and the phase shift between dipole signal, d, 
and laser electric field are plotted. At resonance, light absorption is maximum, and the 
center of charge motion of the electron is approximately 90° out of phase with the laser 
field. 

The resonance in figure Efb) is a collective plasmon resonance. In 3D plasmas and 
clusters the collective excitation of the bound electrons is referred to as a plasma-wave 
and as a plasmon, respectively. We define here the term plasmon as the corresponding 
collective motion of the bound electrons of our ID model system. The collective 
frequency depends on the system geometry, explaining the difference between plasma 
and plasmon frequency. As our model system is neither the ID limiting case of a bulk 
nor of a sphere, the usual 3D plasma/plasmon frequency is not applicable. The ID 
plasmon frequency is determined by the ID geometry of our model system. As an 
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analytical expression is currently not known, we use the above determined numerical 
value. 

Single electron excitations cannot explain the resonance in figure 0] First because 
they have a narrow linewidth. The broad width (« 0.1 at.u.) of the resonance is a strong 
indicator of a collective process. Second, the oscillation shown in figure Eta) decays. The 
decay is also a typical signature of a plasmon, as due to microscopic collisions energy is 
transferred from the collective electron motion into thermal, single electron motion. In 
contrast to that, the life time of a single-electron excitation is infinite. We have tested 
that the decay is not an artifact of the MCTDHF formalism. Increasing the number 
of determinants does not change the time-dependence of the dipole signal significantly. 
Also, this decay does not come from ionization and a decrease of the norm at the 
absorbing boundaries. We find that the ionization yield is virtually constant (increasing 
by 0.003 during the last 40 optical cycles of the simulation time), while the amplitude 
of the plasmon oscillation is reduced by almost a factor of two. 

The match between the plasmon frequency of uj p = 0.35 at.u. and the frequency 
of the dipole oscillation in figure Ufa) proves that the plasmon excitation is responsible 
for the second cut off observed in HHG. Moreover, further analysis shows that multi- 
plasmon states are responsible for the second plateau in the high harmonic signal 
(figure |2J). The oscillation in figure Ufa) has non-sinusoidal components. A Fourier 
transform of d(t) shows frequency components at multiples of the plasmon frequency, 
kuj p , k = 1,2,3,.... Hence, the non-sinusoidal behavior arises from the interference 
between these multi-plasmon states, quivering at multiples of the plasma frequency. 
Consequently, the width of the second harmonic plateau is determined by the highest 
order of the excited multi-plasmon state, which is k = 4 in figure 13 The multi- 
plasmon excitation comes from the nonlinear (anharmonic) part of the binding potential. 
Whereas, close to the center the potential has a quadratic (harmonic oscillator) space 
dependence supporting a single plasmon, the far-range Coulomb part of the potential 
adds nonlinear terms responsible for the creation of multiple harmonics of the plasmon 
oscillation. 

In contrast to single electron systems, in ME systems the collective excitation energy 
adds to the harmonic cut off. The reason is that the collective energy stays in the 
remaining bound electrons and does not get lost while the valence electron makes its 
excursion into the continuum. To elucidate this point we define the ground state energies 
of the neutral and singly ionized system, and Eq 3 \ respectively. The energies of 
the according plasmon states are given by Eq + uj p and E^ +u p . Here, we neglect the 
difference in the plasmon frequencies between the neutral and the singly ionized state, 
since the difference is of the order of the difference between two adjacent harmonics. As 
a result, in both cases the HOMO potential is given by J p = E^ - Ef\ Although for 
our ME system the HOMO ionization potential for the plasmon state is slightly smaller, 
this is a reasonable approximation. In particular, since the difference further decreases 
for an increasing number of electrons and will eventually disappear in real ME system 
which usually have considerably more than four electrons. 
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Figure 5. Energy level diagram for different pathways for HHG in polarizable 
molecules, (a) standard HHG: an electron is set free by tunnel ionization (broken 
arrow), quivers in the electric laser field (thin full arrow) gaining kinetic energy, E\^ lril 
and recombines to the 4-electron ground state, (dashed dotted arrow), (b) as 

(a) but starting from and returning to the plasmon state E^ + lu p . (c) as (b) but 
returning to the ground state, Eq 4 \ Here, Eq 3 ^ and E^ + u> p denote the singly ionized 
ground state and singly ionized plasmon state, respectively. 



There are different pathways by which HHG can take place, which are illustrated in 
figure El Before ionization the system is in its ground state, Eq 4 \ remains in the ground 
state, Eq , after ionization of the valence electron and returns upon recombination to 
its initial state [figure 03(a)]. This is the standard HHG situation. The system may 
also start out in a plasmon state, E^ + u p} remains in the plasmon state, Eq 3 ^ + u p} 
after ionization, and returns to its four electron plasmon state upon recombination 
[figure Efb)]. For both cases the cut off law is 3.17U P + [Eq^ — Eq 3 ^] = 3.17£/ p + I p since 
for the latterthe plasmon frequency chancels out. However, if ionization starts from 
the plasmon state, Eq^ + u p , but the electron returns to the ground sate, E^, upon 
recombination, the plasmon energy is converted into harmonic photon energy, extending 
the cut off, i.e. 3.17U P + I p + u p [figure a)]. 

So far we have discussed the single system response. At the moment, the observation 
of HHG from excited states is experimentally untested. A significant macroscopic signal 
will only be created, when the individual systems emit harmonic radiation in phase. 
Hence, the question has to be answered whether the second plateau can be detected 
in a macroscopic medium consisting of many individual systems. We believe that this 
is the case for the following reason. The phase of the harmonic signal emitted by a 
single system is determined by the phase difference between ground state and excited 
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(plasmon) state. In order for coherent emission to occur, the phase difference between 
ground and excited state has to be the same for all emitters. As the plasmon excitation 
is driven by the laser field, the phase difference between ground and (plasmon) excited 
state is exclusively determined by the laser field, and therefore is the same for each 
system. As a result, the multi-plasmon states are (laser) phase locked to the ground 
state. The contributions from individual atoms add up coherently and HHG can take 
place from the ground as well as from excited states. 

Finally, with respect to experimental observation, we believe that metal clusters 
with their simple geometry and with plasmon energies of a couple of eV are good 
candidates. Although atoms also support collective oscillations, which are known as 
giant (shape) resonances [201; their life time is usually very short, 3-4 times the plasmon 
oscillation period. As a result the plasmon will likely decay before the active electron 
can return and create harmonic radiation. 

3. Conclusion 

High harmonic generation (HHG) in complex multi-electron (ME) systems was 
investigated within the framework of the multi-configuration time-dependent Hartree- 
Fock (MCTDHF) method. Our analysis of HHG spectra in complex multi-electron 
systems revealed two plateaus. The first cut off agrees with the cut off law of noble 
gases. The second plateau is generated due to non-linear excitation of collective plasmon 
oscillations. It arises from electrons that are ionized from an excited plasmon state 
and recombine to the ground state. The ratio of first to second plateau determines 
the population of the plasmon-states. From the extension of the second plateau the 
order of the highest excited plasmon state can be inferred. Thus, the here identified 
plasmon signature presents a novel tool for the investigation of the non-perturbative 
multi-electron dynamics in complex materials, a regime that is experimentally very 
difficult to access otherwise. 
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